## commented out for replication files 

#  install.packages("xtable")
#  install.packages("ggplot2")
#  install.packages("dplyr")
#  install.packages("rmeta")
#  install.packages("rdrobust")

library(xtable)
library(ggplot2)
library(dplyr)
library(rmeta)
library(rdrobust)

## Define treatment dummy and assign weights 
load("het_data.rdata")
data_both <- 
  data_list[[7]] %>%
  mutate(treated = two.days > 0)

## Loop over bandwidths

plot_data <- 
  data_frame(year = rep(c(2009, 2013, 2014, 2015, "Precision weighted \naverage"), each = 6),
             category = rep(rep(c("All", 
                                  "Parent does not \nlive with child", 
                                  "Parent does \nlive with child"), each = 2), 5), 
             type  = rep(c("Conventional RD", "Robust RD"), 15),
             eff   = rep(NA, 30),
             se    = rep(NA, 30),
             place = rep(seq(from = 8, to = 1, -1)[-c(3,6)], 5))

data_close <- 
  data_both %>%
  mutate(voted     = round(par_vote*n),
         abstained = round((1-par_vote)*n)) 

data_voted <- 
  as.data.frame(lapply(data_close, 
                       function(x,p) rep(x,p), 
                       data_close[["voted"]])) %>%
  mutate(voted = 1)

data_abstained <- 
  as.data.frame(lapply(data_close, 
                       function(x,p) rep(x,p), 
                       data_close[["abstained"]])) %>%
  mutate(voted = 0)

data_both <- 
  rbind(data_voted, data_abstained)

years <- c(2009, 2013, 2014, 2015)

for (i in 1:4){
  year_vec <- years[i]
  data_year <-
    data_both %>%
    filter(year == year_vec)
  
  rd <- with(data_year, rdrobust(y = voted, x = two.days))
  
  plot_data[1 + (i-1) * 6, 4:5] <- c(rd$coef[1], rd$se[1])
  plot_data[2 + (i-1) * 6, 4:5] <- c(rd$coef[3], rd$se[3])
  
    for(j in 1:2){
    data <- 
      data_year %>%
      filter(par_child_cohabit == (j-1)) 
    
    rd <- with(data, rdrobust(y = voted, x = two.days))

    plot_data[1 + 2 * j + (i-1) * 6, 4:5] <- c(rd$coef[1], rd$se[1])
    plot_data[2 + 2 * j + (i-1) * 6, 4:5] <- c(rd$coef[3], rd$se[3])
  
    }
  
}

for(i in 1:6){
  subset <- seq(from = i, to = 18 +  i, by = 6)
  plot_data[24 + i, 4] <- meta.summaries(d  = unlist(plot_data[subset, 4]), 
                                         se = unlist(plot_data[subset, 5]))$summary  
  subset <- seq(from = i, to = 18 +  i, by = 6)
  plot_data[24 + i, 5] <- meta.summaries(d  = unlist(plot_data[subset, 4]), 
                                         se = unlist(plot_data[subset, 5]))$se.summary  
  
}

plot_cohab <- 
  ggplot(plot_data, 
         aes(x = eff,
             xmax = eff + qnorm(0.975) * se,
             xmin = eff + qnorm(0.025) * se,
             y = place,
             alpha = type,
             shape = type)) + 
  geom_errorbarh(height = 0) +
  geom_point(size = 2) + 
  facet_grid(. ~ year) +
  theme_classic() + 
  theme(legend.position = "top") +
  scale_x_continuous(name = "") +
  geom_vline(xintercept = 0, linetype = "dashed", alpha = 0.3) + 
  scale_y_continuous("", breaks = c(1.5, 4.5, 7.5), 
                     labels = c("Parent does\nlive with child",
                                "Parent does not\nlive with child", 
                                "All")) +
  guides(alpha = guide_legend(title = ""), shape = guide_legend(title = "")) + 
  scale_alpha_discrete(range = c(0.5, 1))


